The fluid— fluid interface in a model colloid— polymer mixture: Application of grand 
canonical Monte Carlo to asymmetric binary mixtures 
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We present a Monte Carlo method to simulate asymmetric binary mixtures in the grand canonical 
ensemble. The method is used to study the colloid-polymer model of Asakura and Oosawa. We 
determine the phase diagram of the fluid-fluid unmixing transition and the interfacial tension, both 
at high polymer density and close to the critical point. We also present density profiles in the 
two-phase region. The results are compared to predictions of a recent density functional theory. 
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In colloid experiments, hard-sphere-like systems can 
be realized in which the corresponding phase behavior is 
of purely entropic origin. An example is the fluid-fluid 
unmixing transition that is observed in solutions of col- 
loids and non-adsorbing polymers EL H- This transition 
is due to a depletion effect 0, 0, 0,0> @ ■ Each colloidal 
particle is surrounded by a depletion zone from which 
polymers are excluded. When two colloids are close to- 
gether, their depletion zones may overlap, thereby in- 
creasing the free volume available to the polymers. This 
results in an anisotropic pressure exerted by the polymers 
onto the colloids which gives rise to an effective attraction 
between the colloids. Such depletion forces also occur in 
mixtures of large and small hard spheres, but in this case 
it is still debated whether phase separation occurs @ . 

A simple model for colloid-polymer mixtures was first 
introduced by Asakura and Oosawa and later indepen- 
dently by Vrij 0] • In this model (the so-called AO model) 
colloids and polymers are treated as spheres with respec- 
tive radii R c and R p . Hard sphere interactions are as- 
sumed between colloid-colloid (cc) and colloid-polymer 
(cp) pairs, while polymer-polymer (pp) pairs can inter- 
penetrate freely. This yields the following pair potentials 
for the AO model: 



u cc (r) 
u cp (r) 



oo for r < 2R C 

otherwise, 

oo for r < R c - 

otherwise, 



(1) 



= o, 



where r is the distance between two particles. The poly- 
mers thus represent ideal polymer coils with a radius of 
gyration R p which can be realized experimentally in a 9 
solvent. 

The AO model has been the subject of many studies in 
the framework of density functional theories (DFT) 0, 0, 
9]. In particular, these studies yielded the phase diagram 
for a wide range of colloid to polymer size ratios. More- 
over, in the case of the fluid-fluid unmixing transition, 
they predicted interfacial tensions that are roughly one 
thousand times lower than those for simple liquids, in 



agreement with experiments [ToL Hl| . A drawback of the 
latter DFTs is that they are mean field theories and thus 
cannot recover the 3D Ising critical behavior observed, 
for instance, in a recent experiment on a real colloid- 
polymer mixture |ll| . 

Monte Carlo simulations are also well suited to study 
the phase behavior of colloid-polymer mixtures. Recent 
simulations in the Gibbs ensemble were performed to de- 
termine phase diagrams of the AO model ^3] an d also 
of a model that considers non-zero interactions between 
the polymers ^jjj. In the latter study even quantitative 
agreement with experiments was obtained. However, in- 
terfaces are absent in the Gibbs ensemble ^4(, so these 
simulations do not enable investigations close to the crit- 
ical point, nor investigations of the interface in the two- 
phase region. 

The general problem in the simulation of asymmet- 
ric binary mixtures (such as the AO model) is that, by 
displacing or inserting a large particle, overlap with a 
number of small particles will likely result. Such over- 
laps are generally unfavorable and will be rejected in the 
majority of cases. If no special steps are taken, one ends 
up displacing mainly small particles while the large par- 
ticles remain essentially frozen. For asymmetric binary 
mixtures in the canonical ensemble a number of special- 
ized algorithms have been developed that circumvent this 
problem 0, 0] . Unfortunately, it is difficult to obtain 
the surface tension accurately in the canonical ensem- 
ble because it must be derived from the rather small 
anisotropy of the pressure tensor in that case: long wave- 
length interfacial fluctuations are hard to equilibrate in 
the canonical ensemble [l7| . 

In this letter we present a grand canonical Monte Carlo 
method that enables direct simulations of asymmetric bi- 
nary mixtures. We use the method to study the fluid- 
fluid unmixing transition in the AO model. Our method 
consists of collective Monte Carlo moves in conjunction 
with an umbrella sampling technique that was recently 
developed by Virnau and Miillcr 18] . This way, we are 
able to calculate the phase diagram of the AO model close 
to the critical point. At the same time we can calculate 
the interfacial tension 7 in the two-phase region using 
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FIG. 1: Schematic picture of the grand canonical Monte 
Carlo moves used in our simulations. Particles to be inserted 
or removed are shaded grey. See text for details. 

the method of Binder We present a comparison of 
our data to recent DFT results 0, 0, and we show 
that the AO model displays 3D Ising critical behavior. 

The grand canonical ensemble requires that the tem- 
perature T, the volume V and the respective chemical 
potentials {/z p , fi c } of the small (polymer) and large (col- 
loid) particles are fixed. In a standard grand canoni- 
cal move one tries to insert or remove a particle using 
Metropolis sampling [2l|. This approach is not efficient 
for asymmetric binary mixtures because the insertion of 
large particles (colloids) is severely hampered by the pres- 
ence of small particles (polymers). The method that we 
present is aimed to circumvent this problem. The main 
idea is to not transfer particles one-at-a-time, but to 
swap clusters of small particles instead. This stimulates 
the formation of voids, in which a large particle can be 
inserted without producing overlap. 

Figs.n^ and^3 demonstrate the insertion of an addi- 
tional large particle into an asymmetric binary mixture 
currently containing N c large particles. First, a point is 
selected randomly in the mixture and a sphere with ra- 
dius 5 and volume Vg — AttS 3 /3 is drawn around it (see 
Fig. n^.). Let n p denote the number of small particles 
inside the sphere: a particle is inside the sphere when 
the coordinates of its center are inside the sphere. Next, 
we choose a uniform random integer n r from the interval 
< n T < m, with m an integer that will be specified 
later. If n T > n p the move is rejected but if n r < n v , 
n r small particles are randomly selected from the sphere 
out of the n p present (these particles are shaded grey 
in Fig. CIV). The n r selected particles are then removed 
from the mixture and a large particle is inserted at the 
center of the sphere (see Fig. 0)3). The new configuration 
is accepted with probability: 



A + = 



1- 



z c V 



(np)! 



-0AE 



with AE the potential energy difference between the ini- 
tial and the final configuration and {z c , z p } the fugacity 
of large and small particles, respectively. The fugacity z 
is related to the chemical potential via z = exp(/3/i) with 
[3 = l/iksT) and ks the Boltzmann constant. 

The reverse move is illustrated in Figs. [T]p and 
First, a large particle is selected at random and a sphere 
with radius 6 is drawn around the center of this particle. 
Next, a uniform random integer n r from the interval < 
n r < m is chosen followed by the selection of n r random 
locations inside the sphere. These random locations are 
marked as crosses in Fig. Qf!. The selected large particle 
is now removed from the mixture and n r small particles 
are placed on the locations selected before (the newly 
inserted particles are shaded grey in Fig. |Tp). The new 
configuration is accepted with probability: 



n c (n p y.(z p v s y 

z c V (n p + n r )! 



-PAE 



(3) 



N c + l{n v -n t )\ {z p V s ) n * 



(2) 



the notation being the same as in Eq. 

It is straightforward to show that the acceptance prob- 
abilities A + and A- enforce detailed balance, which en- 
sures that the algorithm is not statistically biased [22| . 
The algorithm is also ergodic because single large parti- 
cles have a finite probability of being inserted anywhere 
in the system with one move. Similarly, a small particle 
can be inserted anywhere via a combination of moves: 
for example, by the insertion of a large particle followed 
by the removal of the same large particle. 

In order to apply the method to the AO model, the 
parameters 5 and m still need to be specified. We use 5 = 
R c + R p which is just big enough to contain one colloid 
in a sea of polymers. The integer m must be chosen large 
enough to allow for the formation of voids. In the pure 
polymer phase, the polymer density equals z p because the 
polymers behave like an ideal gas. The insertion sphere 
will then contain z p Vs polymers on average and thus we 
choose m slightly above this value. With this choice of 8, 
the insertion of a colloid can only succeed if all polymers 
are removed from the insertion sphere Vs. This will in 
general happen a fraction 1/m of the time. To boost 
the acceptance rate, we choose to remove all polymers 
from Vs when we attempt to insert a colloid, provided 
their number does not exceed m: moves that attempt to 
remove more than m polymers are rejected. To maintain 
detailed balance, the acceptance probabilities A + and A_ 
must be multiplied by 1/m and m, respectively. 

The phase diagram of the AO model can be expressed 
in terms of the reduced polymer packing fraction 77^ = 
(47r/3)Zp-Rp as a function of the colloid packing frac- 
tion rjc = (47r/3)i?j? N c /V. This is analogous to the 
temperature-density phase diagram for simple fluids. In 
case of the AO model, rj p plays the role of inverse tem- 
perature and rjc that of order parameter. In order to 
determine the coexistence curve of the fluid-fluid unmix- 
ing transition, we calculate for a given value of rf p the 
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FIG. 2: Phase diagram of the AO model with q = 0.8. The 
points are the binodals as obtained from the simulation at the 
indicated volumes of the simulation box. The solid line is the 
binodal from DFT 



probability distribution P(jj c ). This is the probability of 
observing a mixture with colloid packing fraction rj c . If 
phase separation occurs P{r] c ) is bimodal: the peak at 
low rj c corresponds to the colloid vapor phase, the peak 
at high r) c to the colloid liquid phase, and the region in 
between is the phase-separated regime. To ensure phase 
coexistence, the colloid fugacity z r is tuned such that the 
area under both peaks is equal |2l) . 

A crucial point in our simulation is the use of a new bi- 
ased sampling technique called successive umbrella sam- 
pling. This technique was recently developed by Virnau 
and Miiller [l^. Combination of our grand canonical 
Monte Carlo move and successive umbrella sampling en- 
ables us to sample P(f] c ) also in regions where, due to the 
free energy barrier separating both phases, P(f] c ) may be 
very low. 

In the following we consider an AO mixture with a 



polymer to colloid size ratio q = R p /R c = 0. 



Note 



that all lengths are given in units of R c = 1. From pre- 
vious studies we expect fluid-fluid phase separation for 
q = 0.8 in a wide range of rf p . The simulations were 
performed in an elongated box with aspect ratio 1/2 and 
periodic boundary conditions. We have performed sim- 
ulations for two different volumes of the simulation box: 
Vi = 4652.4 and V 2 = 9314.9. When the colloid packing 
fraction reaches 0.45, the smaller volume contains 500 
colloids and the larger volume 1000 colloids. 

Fig-l2shows the phase diagram for the two system sizes 
V\ and V<x. Comparison of the two data sets shows that 
finite-size effects are relatively small, even close to the 
critical point. Also included in Fig. (21 is the phase dia- 
gram obtained from recent DFT [a, |9j . We observe that 
DFT underestimates the location of the critical point 
by about 30%. More importantly, DFT yields the typ- 



FIG. 3: Reduced interfacial tension 7* = 4i?^7 as a function 
of the difference between the colloid packing fractions in the 
coexisting liquid (L) and vapor (V) phases. The open circles 
represent our simulation data. The solid line is a DFT re- 
sult |2Cft . The inset shows the logarithm of the probability 
distribution P(?? c ) for rf v = 1.225 with Fl defined in the text. 



ical mean-field parabolic shape of the binodal (critical 
exponent /3 = 1/2) while the simulation yields the ex- 
pected flatter binodal (/3 « 0.325; Ising model univer- 
sality class [23j). We will demonstrate below that the 
simulation indeed displays Ising critical behavior. 

The probability distribution P(r) c ) can also be used to 
extract the interfacial tension 7 between the coexisting 
colloid vapor and colloid liquid phases. To this end one 
uses a formula that was first derived by Binder |l9| : 



L — »oo 



P max (r? c ) 
P min (Vc) 



(4) 



with P max {r) c ) and P min (r] c ) the value of P(i lc ) at its 
maxima and its minimum, respectively, and L the length 
of the simulation box parallel to the interface. An ex- 
ample distribution is shown in the inset of Fig. for 
r/p = 1.225. Note that the presence of a flat region be- 
tween the two peaks is important for an accurate esti- 
mate of 7. This is enforced in the simulation by using an 
elongated box. 

Fig.[3]shows the reduced interfacial tension 7* = 4R 2 j 
as a function of the difference between the colloid packing 
fractions in the coexisting phases, as obtained from our 
simulation together with a recent DFT result. Again, the 
DFT result deviates from the simulation by about 30%. 
This demonstrates that the rather perfect agreement of 
DFT with experimental data as claimed in Ref. p] is co- 
incidental. According to our simulation data, the values 
of 7* for the AO model underestimate experimental data 
significantly which shows that more sophisticated mod- 
els are required to describe colloid-polymer mixtures on 
a quantitative level. 
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FIG. 4: The width of the colloid interface W/ (2R C ) as a func- 
tion of the relative distance from the critical reduced polymer 
packing fraction for two different system sizes as indicated. 
The solid lines are fits to the 3D Ising power law described 
in the text. The inset shows three actual colloid density pro- 
files along the direction perpendicular to the interfaces for 
»7p = {0.82,0.98, 1.20}. 



We have carried out additional simulations to calculate 
colloid density profiles in the two-phase region. These 
simulations were performed using 1156 and 2889 colloidal 
particles at a colloid packing fraction of 0.13. Periodic 
boundary conditions were again used and the aspect ra- 
tio of both boxes was 1/3. This corresponds to a box 
length of L z = 69.46 and L z = 94.27 for the smaller and 
larger system, respectively (here L z is the box length 
perpendicular to the interface). Three density profiles 
for the larger system are shown in the inset of Fig. 0] We 
have estimated the (10%-90%)-interfacial width W from 
these profiles by fitting the profiles to a hyperbolic tan- 
gent. Since the interfacial width is expected to diverge 
with the same exponent as the bulk correlation length 
near the critical point, one expects W oc (rf^ — rf p crit ) 
with v — 0.63 corresponding to 3D Ising critical behav- 
ior [2i|. Fig. Q] shows that the data for both system sizes 
is consistent with this power law. We can also infer from 
Fig. |1| that W depends strongly on L z , even far away 
from the critical point. This is most likely due to capil- 
lary waves [24J and care has to be taken when comparing 
W from this simulation to interfacial widths obtained in 
experiments or analytical theories. 

In summary, we have presented a grand canonical 
Monte Carlo method which is well suited to simulate 
asymmetric binary mixtures. It is particularly powerful 
when combined with a re-weighting scheme: both the 
phase diagram and the surface tension can be obtained 
accurately in that case. We have used the method to de- 
termine the coexistence line of the fluid-fluid transition 
in the AO model with high accuracy. We have also pre- 



sented new analysis of the interface between coexisting 
phases in the AO model, in particular estimates of the 
interfacial tension. 
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